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ABSTRACT 


The  present  report  describes  research  accomplished  during  the 
first  year  of  the  cooperative  program  aimed  at  determining  the  nature 
of  the  seismic  radiation  from  small  strike-slip  earthquakes.  Our 
activities  have  been  in  two  areas,  to  develop  and  operate  four  long 
period  portable  stations  in  the  area  to  be  monitored,  and  to  develop 
the  necessary  ‘-.heoretical  framework  for  the  interpretation  of  the  ob¬ 
servational  data. 

In  summary,  during  the  first  year  contract  period  we:  (1)  de¬ 
veloped  four,  low  power  multi -component  long  period  trailer  units  with 
broadband  recording  capability  employing  both  analog  and  digital  record¬ 
ing  capability  (7  channels  each);  (2)  upgraded  eight  existing  trailer 
units  to  broadband,  analog  recording  capability  to  record  both  long  and 
short  period  data;  (3)  installed  and  operated  a  four  trailer  array  in 
Bear  Valley,  California  area  and  successfully  recorded  two  small  events 
(magnitudes  4  and  3.5)  from  Bear  Valley  and  numerous  teleseismic  events; 
(4)  investigated  near  and  far  field  radiation  from  relaxation  source 
models  of  earthquakes  and  extended  the  theory  to  a  variety  of  source 
geometries  and  prestress  conditions.  It  was  found  that  the  low  frequency 
behavior  of  the  near  fi el d  was  1/f,  as  expected.  The  low  frequency  be¬ 
havior  of  the  far  field  varies  between  a  flat  spectrum  and  a  spectrum 
decreasing  as  f  with  decreasing  frequency  away  from  a  spectral  peak, 
depending  on  whether  the  prestress  field  is  uniform  to  infinity  or  con¬ 
centrated  in  a  zone  of  characteristic  dimension  equal  to  a  few  times  the 
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fault  dimension.  Larger,  yet  finite,  prestress  zones  produce  a  broad 
spectral  peak,  the  spectrin  remaining  nearly  flat  with  decreasing  fre- 
quency  to  a  characteristic  frequency  which  is  controlled  by  the  dimen¬ 
sion  of  the  prestress  zone,  below  which  the  spectrum  begin?  to  occrease 
2 

again  as  f  for  lower  frequencies.  Within  a  wavelength  or  so  from 
the  source  the  spectrum  is  dominated  by  the  near  field  spectral  component. 
At  high  frequencies,  above  the  characteristic  "peak"  or  "comer"  fre¬ 
quency  determined  by  the  rupture  length  and  rupture  velocity,  the  spec¬ 
trum  decreases  as  1/3.;  (5)  investigated  near  field  wave  propagation 
in  a  layered  medium  using  the  Cagnai rd  method  incorporating  various  simple 
source  models  and  applied  the  method  to  the  prediction  of  the  field  from 
several  small  earthquakes,  wi th  reasonable  first  order  agreement  with 
observations;  (6)  incorporated  complex  source  models  (e.g.,  relaxation 
and  dislocation  types)  in  surface  and  body  wave  (ray  theoretic)  programs. 
Both  near  and  far  field  source  terms  are  included  and  these  computations 
should  be  accurate  in  the  teleseismic  and  intermediate  distance  ranges 
(up  to  a  wavelength  or  two  from  the  source).  The  programs  are  being  used 
to  predict  and  from  earthquakes .  Preliminary  results  for  m^ 
vs.  Mg  are  in  overall  agreement  with  observations.  A  cutoff  m^  value 
is  predicted  implying  that  if  1  cps  energy  is  used  to  calculate  m^  , 
then  the  maximum  for  any  earthquake,  however  large,  will  be  about 

7.0,  and  for  "normal"  prestress  levels  of  a  few  hundred  bars,  the  maximum 
will  be  near  6.0.  An  Mg  cutoff  value  of  around  10  to  11  is  also 
predicted. 
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I.  Introduction 

The  purpose  of  the  research  conducted  under  this  contract  is  to 
determine  the  detailed  nature  of  the  seismic  radiation  spectrum  from 
small  earthquakes,  especially  the  nature  of  the  long  period  part  of  the 
spectrum.  The  objective  is  to  not  only  verify  discrimination  criteria 
for  various  types  of  small  earthquakes  (e.g.,  mb  vs.  criteria)  but 
to  obtain  a  fundamental  understanding  of  why  such  criteria  work.  This 
implies  that  we  must  obtain  a  model  for  an  earthquake  that  is  sophisti¬ 
cated,  incorporating  the  basic  physics  of  the  process,  and,  in  addition, 
that  we  must  be  able  to  predict  the  near  field  wave  propagation  effects 
as  well  as  far  field  wave  propagation  in  order  to  interpret  the  observa¬ 
tions  in  terms  of  the  earthquake  model. 

Our  work  under  this  contract  during  the  first  year  has  therefore 
been  divided  into:  (1)  Field  operations  and  instrumental  systems  fabri¬ 
cation.  (2)  Investigations  of  wave  propagation  theory  from  complex 
seismic  sources  in  the  near  and  far  field  distance  ranges.  And  (3)  Model 
ing  of  earthquake  by  numerical  and  analytical  relaxation  sources. 

In  addition,  of  course,  data  reduction  and  interpretation  would  be 
an  essential  part  of  our  program;  however  no  well  recorded  earthquake 
was  obtained  during  the  present  contract  period,  so  that  no  extensive 
interpretive  work  was  done.  Hopefully  we  will  have  the  opportunity  to 
analyze  an  earthquake  from  the  Bear  Valley  area  during  the  course  of  the 
next  contract  period. 

In  the  following  sections  we  discuss  our  accomplishments  in  some 
detail,  with  many  of  the  tharonetical  results  incorporated  in  this  report 
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prior  to  publication  in  order  that  other  program  investigators  be  aware 
of  the  implications  of  our  work  in  their  own  investigations. 

1L  Fie1d  Program  and  Instrumentation 

Our  part  of  the  cooperative  field  program  was  to  provide  record¬ 
ing  of  the  long  period  radiation  from  small  earthquakes  in  the  magnitude 
range  3.0  to  5.0.  Four  systems  were  to  be  designed  and  placed  in  the 
designated  field  site  near  Bear  Valley,  California.  In  this  section  we 
describe  the  system  design  used  to  achieve  broad  band  recording  in  the 
range  from  a  few  cycles  out  to  periods  of  around  60  seconds,  as  well  as 

the  field  program  actually  undertaken  in  cooperation  with  other  investi¬ 
gators. 

0)  Instrumental  Characteristics  of  Mark  I  and  Mark  II  Trailer 
Uni  *  s 

In  order  to  achieve  a  highly  portable  and  reliable  field  record¬ 
ing  system  for  long  period  seismic  radiation  from  small  events,  we  chose 
to  employ  internally  recording  trailer  units,  using  a  relatively  short 
period  seismometer  (with  adjustable  period  from  5  to  15  sec)  in  order  to 
avoid  drift  problems  associated  with  longer  period  seismometers.  Further¬ 
more,  use  of  the  5  sec  seismometer  minimizes  parasitic  effects  inherent 
in  long  period  recording,  wherein  high  frequency  pulses  result  in  non¬ 
linear  response  producing  long  period  motion  of  the  recorder.  In  order 
to  extend  the  useful  response  of  the  system  to  long  periods,  we 
apply  a  variety  of  high  frequency  filtering  followed  by  amplification  of 
the  low  frequency  output  using  very  low  noise  amplifiers. 
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Two  separata  trailer  systems  are  employed  in  the  present  program, 
both  of  then:  are  self  contained  recording  systems  with  their  own  battery 
and/or  solar  power  available  if  necessary.  The  newest  (Mark  II)  units, 
constructed  under  the  present  program,  record  both  digitally  and  in 
analog  form  with  a  rather  wide  variety  of  filter  settings 
available,  resulting  in  both  long  period  and  short  period  recording  capa¬ 
bility.  The  digital  recording  is  used  for  the  long  period  data,  and  with 
a  sample  rate  of  2  samples/sec  will  record  seven  channels  of  data  for 
approximately  8  days  before  a  tape  replacement  is  necessary.  Two  differ¬ 
ent  filter  settings  can  be  used,  or  two  different  gain  levels  can  be 
recorded.  The  analog  recording  can  be  carried  on  simultaneously  and  can 
be  used  for  either  high  or  low  frequency  recording  (using  variable  fil¬ 
tering)  of  the  5  sec  seismometer  output  or  it  can  be  used  to  simultaneously 
record  data  from  a  different  type  of  seismometer.  The  analog  system 
records  over  a  40  db  range. 

Since  the  construction  of  the  new  Mark  II  units  was  in  progress 
over  much  of  the  current  report  period,  we  employed  four  of  our  older 
(Mark  I)  trailer  units  as  recorders  in  the  field  over  this  period.  We 
have  eight  such  units,  all  of  which  record  on  FM  magnetic  tape  and  on 
film.  We  modified  all  eight  units  to  function  as  broad  band  systems,  and 
in  particular,  when  using  the  5  sec  seismometers,  the  system  response  is 
very  close  to  the  long  period  response  of  the  Mark  II  units.  These  units 
will  be  replaced  in  the  field  by  the  new  Mark  II  systems  as  they  become 
available.  At  this  time  the  Mark  II  units  are  complete  and  are  being 
field  tested,  so  that  replacement  should  begin  soon.  (One  or  two  of  the 
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Mark  I  units  may,  however,  be  left  in  the  field  to  provide  additional 
coverage  of  the  monitored  site  provided  we  can  afford  the  costs.) 

The  system  response  of  the  Mark  I  and  II  portable  trailer  units 
have  been  obtained  by  shake-table  tests  and  are  given  in  Figures  1 
through  4.  Figure  1  shows  the  pass  band  and  gains  available  from  the 
amplifiers  used  in  these  systems.  Figure  2  gives  the  amplitude  and 
phase  response  for  the  Mark  I  systems  operating  either  as  broadband 
recorders  or  as  "long  period"  recording  systems.  To  date  we  have 
found  that  the  Mark  I  systems  function  exceptionally  well  at!  long  period 
recording  systems.  In  the  next  section  (1 1-2)  we  show  two  small  earth¬ 
quakes  (m^  -  4  and  m^  -  3.5)  from  the  Bear  Valley  area  as  examples 
of  this  recording  capability. 

Figure  3  shows  the  system  response  of  the  Mark  II  units  with 
third  order  low  pass  filtering  at  the  various  filter  settings  available. 
Curve  C,  with  the  filter  setting  at  20  sec  is  equivalent  to  Curve  B  in 
Figure  2  which  shows  the  long  period  response  of  the  Mark  I  units.  This 
response  is  what  we  are  currently  using  in  the  field. 

Figure  4  shows  the  response  of  two  Mark  II  systems  using  fourth 
order  low-pass  filtering  and  the  available  variety  of  pass  bands.  Either 
the  third  or  fourth  order  filters  can  be  used. 

(2)  Station  Configuration  and  Events  Recorded 

Figure  5  shows  the  configuration  of  the  Bear  Valley  array  along 
with  a  general  identification  of  instruments  employed  and  field  program 
participants.  This  array  has  been  completely  operational  for  only  the 
last  4-5  months  and  in  that  period  no  events  occurred  within  the  central 
monitored  area,  nor  in  fact  within  the  general  area  enclosed  by  the  most 


Figure  3 
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extremely  separated  elements  of  the  array.  However,  prior  to  the  com¬ 
pletion  of  the  array  two  small  events  occurred  which  were  recorded  by 
our  long  period  trailers.  Even  though  we  cannot  analyze  these  events 
using  the  full  array,  the  recordings  provide  a  demonstration  of  the 
ability  of  these  portable  units  to  adequately  record  long  period 
radiation  from  events  in  the  magnitude  range  of  interest.  Figure  6  shows 
the  two  events  as  recorded  by  one  of  the  trailer  units. 

It  is  felt  that  the  units  presently  in  the  field  will  record 
mb  s  3-5  events  out  to  the  most  extreme  distance  and  that  the  Mark  II 
units  will  be  somewhat  more  sensitive. 

III.  Theoretical  Developments 

In  anticipation  of  successfully  recording  earthquakes  in  the  near 
and  far  field  distance  ranges  during  this  program  we  have  developed  the 
theoretical  framework  for  the  interpretation  of  near  and  far  field  wave 
propagation  effects  as  well  as  a  more  comprehensive  relaxation  source 
representation  which  treats  near  field  effects  in  detail. 

The  problems  faced  in  treating  the  near  field  of  an  earthquake  are 
two-fold.  First:  An  earthquake  is  a  volume  source,  in  that  energy  is  re¬ 
leased  from  within  some  finite  volume  resulting  in  a  net  decrease  in 
stored  strain  energy  within  that  region;  the  question  then  is  whether 
measurements  of  the  radiation  field  made  withir  this  source  volume  will 
be  fundamentally  different  than  measurements  made  when  outside  of  it. 

In  the  near  field  it  is  quite  literally  the  case  that  the  real  source  is 
distributed  in  space  all  around  the  point  of  observation  instead  of  being 
localized  at  some  point,  or  on  some  surface.  In  this  case,  the  usual 
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theoretical  treatments  of  energy  radiation,  which  assume  a  localized  source, 
are  inadequate,  or  at  least  have  implicit  in  their  formulation  an  assumption 
regarding  the  spatial  origins  of  the  energy  release  which  are  extremely  restrictive. 

Thus  the  problem  is  to  determine  the  field  to  be  observed  at  a 
point  which  is  within  the  source  itself;  in  the  case  of  an  earthquake, 
within  the  zone  of  stress  relaxation. 

The  second  problem  is  that  conventional  "far  field"  approximations 
used  in  wave  propagation  theory  must  be  given  up  and  the  propagation  of 
energy  in  the  rather  complex  medium  must  be  calculated  more  exactly  for 
near  field  observations.  Furthermore,  interference  effects  related  to 
multiple  reflections,  and  hence  waveguide  phenomena  giving  rise  to 
ordinary  surface  waves,  for  example,  cannot  be  fully  effective  in  the 
near  source  region  since  the  required  constructive  or  destructive  inter¬ 
action  of  reflected,  refracted  and  diffracted  waves  will  not  have  completely 
occurred.  This  means  that,  in  the  very  near  source  region,  it  is 
necessary  to  calculate  the  field  using  a  Cagniard  (or  generalized  ray 
theory)  method  which  takes  explicit  account  of  individual  generalized 
rays  and  sums  them  to  provide  the  predicted  radiation  field.  Alternatively, 
the  so-called  "leaking  mode"  theory  can  be  used,  or  a  full  numerical 
calculation  using  finite  difference  or  finite  element  methods. 

In  addressing  ourselves  to  these  problems,  we  decided  in  the  first 
jlace  that  the  problems  were  such  that  without  special  attention  devoted 
to  them,  one  could  not  hope  to  obtain  a  meaningful  interpretation  of  the 
field  data  to  be  obtained  in  this  program.  Therefore  we  began  Investigations 
in  three  separate  areas;  (1)  on  the  theory  of  the  source  in  the  near  field 
range  for  the  receiver  inside  the  prestress  zone  as  well  as  outside  the 


zone; 
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(2)  on  the  inclusion  of  complex  relaxation  source  models  in  surface 
wave  propagation  programs  wherein  all  near  field  effects  are  included; 

(3)  and  on  the  inclusion  of  complex  source  models  in  the  Cagnaird 
genralized  ray  theory,  again  retaining  the  near  distance  radiation 
effects  from  the  source.  Further,  in  conjunction  with  other  work  related 
to  earthquake-explosion  discrimination,  we  began  systematic  investigations 
of  source  characteristics  utilizing  these  theoretical  capabilities. 

In  the  following  subsections  (II-l  through  II-3)  we  summarize, 
in  some  detail,  the  results  of  the  theoretical  investigations.  In 
Section  IV  we  give  a  summary  of  some  of  the  preliminary  results  obtained 
through  applications  of  the  theory  to  the  discrimination  problem, 
particularly  in  terms  of  source  characteristics  which  lead  to 
discrimination  on  the  basis  of  vs  Mg. 


(!)  Extensions  of  Relaxation  Source  Theory  for  Near  Field  Representations 

The  solution  can  be  expressed  in  terms  of  potentials.  Archambeau  (1972) 
obtained  (eq.  5-1). 


xa  a  =  i...4  represent  the  three  components  of  rotation  and  the  dilatation. 

Va  is  the  appropriate  wave  velocity,  x*  are  the  equilibrium  fields  appropriate 
to  a  moving  boundary  X  in  the  medium.  l^(t  )  is  the  volume  exterior  to  E.  as  a 
function  of  t^,  the  source  time,  and  t  is  the  total  rupture  time. 

For  the  case  of  the  self  similar  problem,  where  £  is  an  expanding  spherical 
rupture,  V(tQ)  is  the  volume  exterior  to  that  sphere. 

However  for  the  case  of  a  propagating  spherical  rupture,  we  *t  .ll  use  an 
approximation,  necessitated  by  the  lack  of  symmetry  of  the  problem.  In  particular, 
we  shall  consider  the  case  of  a  unilateral  rupture  and  consider  a  spherical  rupture 
propagating  from  the  point  O'  along  the  axis,  with  a  radius  R(t  ), 

and  center  at  d(tQ)  on  this  axis  (Figure  1).  Because  of  theoretical  complications 
associated  with  healing  phenomena  (after  rupture),  the  treatment  given  in  this 
report  will  be  confined  to  the  case 

d(to)  <  R(to) 


More  complicated  rupture  geometries,  of  greater  generality,  will  be  treated 

later.  However,  this  enveloping"  model  will  show  all  the  essential  characteristics 
of  importance. 
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In  order  to  preserve  the  symmetry  of  the  problem,  the  volume  of 
integration  l/(to)  is  then  defined  as  the  volume  exterior  to  the  sphere 
8(tQ)  (Figure  1),  centered  at  the  hypocenter,  and  going  through  the 
rupture  front.  The  radius  of  this  sphere  is  then  VRtQ.  This  neglects 
contributions  from  inside  8(tQ)  but  the  approximation  is  justified  by  the 
fact  that  all  of  the  non-elastic  phenomena  associated  with  failure  take  place 
within  8(tQ),  and  the  elastic  energy  stored  inside  B(tQ)  can  safely  be  assumed 
to  be  dissipated  by  these  phenomena. 

It  is  convenient  to  take  Fourier  transforms  of  (1)  with  respect  to  t,  and  then 


we  have 


a  Jo  Jo 

v<to) 


*  -ik  r* 
a 

a  e 


dr  dt 


-ikar* 

The  function  e  /r*  is  the  Green's  function  for  the  Helmholtz  equation 

in  the  infinite  domain.  Its  use  in  equation  (2)  is  an  approximation:  The 
solution  given  by  (2)  will  satisfy  the  initial  value  conditions,  but 
Xa(r,t)  will  not  necessarily  satisfy  the  appropriate  boundary  conditions 
on  Z  for  all  t.  In  order  to  satisfy  those  boundary  conditions,  one  would 
have  to  superpose  general  solutions  of  the  homogeneous  equations  and  (2).  This 
additional  term  would  in  essence  represent  the  interaction  of  the  dynamic 


fields  with  the  boundary  E,  in  other  words,  the  scattered  fields.  Ignoring 
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the  scattered  fields  is  equivalent  to  making  the  inclusion  within  E 
transparent  to  incoming  waves  and  can  indeed  be  taken  as  the  definition  of 
transparency.  To  model  underground  explosions,  we  consider  an  expanding 
sphere  of  radius  R(tQ),  upon  which  the  tractions  vanish,  center  at  O'.  The 


initial  values  x*  are  given  by 


■  (r‘>*  5  (‘ 


cos  m  d> 
2m  Tc 


+  b^0^  sin  m  <|>  )  (cos  0  ) 
2m  o/2  o 


Where  the  coefficients  a^“^  and  b^^  are  given  by  Archambeau  (1972,  eq.  4.9 


and  4.10) 
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where  o^j  represents  the  prestress,  chosen  to  be  pure  shear  and 
homogeneous  for  this  case. 

On  the  other  hand,  in  the  case  of  a  propagating  rupture,  one  has  to 

take  into  account  the  fact  that  at  time  t  ,  the  failure  zone  is  centered 

o 

on  the  axis  at  a  distance  d(tQ)  from  o'  .  Using  addition  theorems  for 
the  Legendre  associated  functions  one  gets  in  that  case 

2 

£ 


v*(r  ,t  )  ■  ~  M  (  a^  cos  m  tf>  +  b^  sin  m  a  \ 

*a Vio’  o  r3  m-0  y  2m  vo  2m  yo  J 


(2-nH-s)j  d(to) 

s!(2-ra)»  r 
s-u  L  o 


*2+s  (cos  eo) 


(5) 


where  the  coefficients  a^  (t  )  and  b^a\t  )  are  again  given  by  (4.1)  and 

2m  o  2ra  o 

(4.2).  Note  that  only  the  harmonics  of  degree  2  are  present.  This  is  not 

the  case  for  non-spherical  ruptures. 

Therefore  we  see  that  (3)  is  merely  a  particular  case  of  (5)  where 

the  sum  over  the  index  s  is  reduced  to  its  first  term  s  *  0.  (5)  is  only 

valid  for  r  <  d(t  ),  however,  this  is  always  the  case  for  r  e  ^(t  ). 

o  o  — o  o 

The  initial  value  X*(r  >t  )  vanishes  like  -^r  at  r  -*■  °°.  However, 

a  — o  o  jo 

r 

if  only  because  of  the  finiteness  of  the  earth,  and  because  the  present 
oj^can  hardly  be  homogeneous  at  very  large  distances  from  the  source, 
one  may  assume  that  this  initial  value  becomes  vanishingly  small  beyond 
a  relaxation  radius  R  .  The  simplest  way  to  approximate  such  a  behavior 
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is  to  use  the  initial  value  given  by  (3)  or  (5)  in  (2),  but  to  truncate 
the  volume  integral  in  (2)  at  the  radius  R  .  The  volume  within  R  will 

O  b 

be  henceforth  referred  to  as  the  relaxation  zone.  Two  possible  geometries 

then  arise.  In  a  first  instance,  the  relaxation  zone  may  be  confined 

to  the  vicinity  of  the  failure,  so  that  most  observer's  points  £  would  lie 

outside  of  it.  This  is  the  geometry  investigated  previously  by  Archambeau 

(1968).  Archambeau  showed  that  the  error  on  the  energy  released  is 

negligible  provided  that  Rg  is  equal  to  about  five  times  the  fault  length 

dQ.  This  case  we  shall  consider  as  one  possible  extreme  behavior  -  the 

other  geometry  corresponds  to  the  case  where  R  may  be  large  and,  in  the 

b 

limit,  infinite  -  observer  points  then  lie  within  the  relaxation  zone. 

This  constitutes  another  extreme  behavior. 

These  two  extreme  behaviors  will  allow  us  to  place  bounds  on  the 
spectral  content  of  the  radiation  fields,  the  reality  being,  of  course, 
between  these  extremes. 

We  shall  now  treat  the  case  of  a  propagating  rupture,  with  a  finite 
relaxation  zone  and  an  observer  inside  this  relaxation  zone.  This  is  the 
most  general  case,  and  the  results  appropriate  to  the  extreme  geometries, 
as  well  as  the  results  for  non-propagating  ruptures  may  be  deduced  easily 
thereafter. 

We  first  rewrite  equation  (2)  in  the  form: 
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Equation  (6)  is  obviously  valid  for  any  of  the  cases  considered  above.  We 
shall  make  use  of  the  usual  spherical  wave  expansions 
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where  the  upper  pair  of  Bessel  functions  are  to  be  used  when  |r|  >  |r  |  , 
and  the  lower  pair  when  [r^J  >  |_r|  > 

Here 
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where  6^  is  the  usual  Kronecker  delta.  We  shall  also  use  the  integral 
relation 
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// 
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Using  (7),  (8),  (9)  in  (6),  inserting  (5)  for  xq  where  appropriate,  and 


splitting  the  last  integral  in  (6)  as 


-  R  „  r  .  R 
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where  the  multipole  coefficients  ^ A 2m»  £C2m,  £D  ,  are  given  by: 
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For  the  non-propagating  rupture,  because  symmetry  is  preserved,  only  the 
term  £  ■  2  survives  in  (10),  (11)  and  (12). 

Despite  its  symmetry,  solution  (10)  is  not  particularly  convenient 
since  the  multipole  coefficients  depend  on  r,  the  hypocentral  distance 
of  the  observer.  A  more  convenient  form  for  computation  can  be  obtained 
by  evaluating  the  integrals  in  (11)  and  (12). 
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It  can  be  shown  (Minster,  1973)  that 
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Inserting  these  formulae  into  (11)  and  (12)  one  gets 
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One  sees  immediately  that  the  only  integrals  left  to  evaluate 


^  0  O  L 


3(t  )  d£"2(t  )  dt 
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which  can  be  evaluated  in  closed  form  (Minster  and  Archambeau,  1973)  and 


J°' 

0 

for  which  closed  forms  can  be  found,  but  which  is  nevertheless  best 
computed  numerically  as  a  finite  Fourier  transform. 

We  can  now  combine  the  r  dependent  terms  in  (13)  and  (14),  and  make 
use  in  (10)  of  the  Wronskian  relation 
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This  allows  us  to  write  the  solution  in  the  final  form. 
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where 
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The  analysis  proceeds  in  an  identical  fashion  for  the  non-propagating 
case,  by  keeping  only  the  £  -  2  terra.  The  second  term  in  equation  (15)  is 
a  non  propagating  term;  in  fact,  it  represents  the  creation  of  the  initial 
value.  To  see  this  take  the  Fourier  transform  of  (5) 
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which  by  simple  redefinition  of  the  indices  can  be  rewritten 

^(“U)  =  S  m?0  P”(a°S  9)  ((iTTT^1)  (  C°S  ^  +  ^  Sin  m4>)  U7) 

because  (15.)  gives  us  the  relative  field,  (measured  with  respect  to 
the  final  equilibrium),  one  has  to  subtract  (17)  from  (15)  to  obtain 
the  observed  radiation  field. 

We  are  now  in  a  position  to  investigate  the  extreme  cases  presented 
above . 

1)  If  we  take  Rg  infinite,  then  we  have  immediately 
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•  P™(cos  0) 

2)  In  the  case  of  a  limited  relaxation  zone  R  <  r  we  need  only  keep 

O 

the  first  term  in  (10)  and  change  the  upper  limit  of  the  integral  over 

rQ  in  (11)  to  Rg  -  consequently  the  solution  has  the  same  form  as  (18) 

(2) 

and  the  integral  1^  (oo)  appearing  in  (16)  has  to  be  replaced  by 
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This  is  specifically  the  case  investigated  by  Archambeau  (1968). 
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furthermore  we  can  write 
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I.  Farfield  approximation  k^r  <<  1 

In  this  case  we  approximate  the  Hankel  function  by  its  far  field  term 
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The  spectral  shape  is  then  easily  obtained,  because  only  the  quadrupole  term 
is  important  for  long  periods. 
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for  to  <<  1 


This  result  is  the  one  obtained  by  Archambeau  in  previous  studies.  It  implies 


that  the  displacement  spectrum  exhibits  at  least  one  peak.  The  displacement 


spectral  density  vanishes  as  u)  at  long  period. 


2)  Rg  =  °° 
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This  gives  a  spectral  shape  similar  to  that  suggested  by  various  authors 
(e.g.  Aki,  1970  and  Brune,  1970),  with  a  "flat"  long  period  level. 


II.  Near  field  behavior 


In  the  near  field,  the  most  important  term  in  the  Hankel  function  is 
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In  this  case  however,  the  displacement  spectrum  is  worth  a  more  detailed 
investigation.  We  have  in  general 


u(r,to)  =  -  V0(r,to)  +  V  xg(r,to) 


where  the  first  term  represents  the  P-wave  radiation,  and  the  second  the 
S-wave  radiation.  As  an  example  we  shall  consider  in  detail,  the  radial 
component  of  motion. 
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We  see  that 
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Such  asymptotic  behavior  would  seem  to  indicate  that  both  P  and  S  waves 
carry  infinite  energy  at  very  long  periods.  Of  course,  the  reason  for 
this  surprising  behavior  is  that  one  cannot  define  P-  and  S-waves  in  the 
static  limit,  and  that  the  mathematical  separation  indicated  above  is  not 
physically  realized,  that  is  P  and  S  energy  overlaps  and  adds  in  the  long 
period  limit.  In  particular,  we  can  show  analytically  that 

ur0r,u))  =  u^p^  +  u<S>  =  for  0)  «  1 

The  proof  is  rather  tedious  but  tractable.  More  importantly,  numerical  results 
exhibit  the  same  behavior  to  a  high  degree  of  accuracy,  and  this’ provides  a 
useful  check  on  the  stability  and  accuracy  of  numerical  calculations. 


High  frequency  behavior 


One  can  show  that  the  asymptotic  behavior  of  the  displacement  spectra 
for  all  the  cases  considered  in  this  report  is  u^)  =  to  +  ». 

Numerical  calculations  also  show  this  asymptotic  behavior.  However 
the  preliminary  investigations  show  that  the  spectra  are  quite  complicated 
at  high  frequency,  due  to  the  significant  excitation  of  higher  order 
multipoles.  The  same  complexity  is  present  in  theoretical  radiation  patterns 
as  well.  Because  analytical  investigations  become  extremely  cumbersome 
for  the  intermediate  frequency  range,  one  has  to  rely  on  numerical  parameter 
studies  to  comprehend  the  nature  of  the  seismic  radiation  in  that  range.  Such 
studies  are  currently  underway. 


Summary 


We  solved  the  problem  of  the  relaxation  of  a  prestressed  elastic  medium 
due  to  a  propagating  rupture  for  two  extreme  prestress  conditions.  The  first 
one,  where  the  relaxation  zone  is  limited  in  size,  may  be  used  to  model  the  effect 


of  a  local  concentration  of  stress.  The  other  one,  where  an  infinite  medium 
is  allowed  to  relax,  models  cases  of  high  stress  levels  on  a  regional  scale, 
and  is  particularly  useful  in  the  investigation  of  near-field  effects.  On 
the  basis  of  these  two  extremes,  one  can  bracket  the  long  period  behavior 
of  the  far  field  displacement  spectrum 

0(w2)  <  [u]  <  0(1)  for  a)  «  1 

The  lower  bound  corresponds  to  the  source  model  previously  proposed  by 
Archambeau  (1968),  and  yields  a  peaked  spectrum.  The  upper  bound  gives  a 
long  period  spectral  shape  similar  to  that  proposed  by  Aki  (1970)  and 
Brune  (1970) . 

The  "near  field"  is  defined  as  that  part  of  the  radiation  field  decaying 

with  distance  as  — —  ,  n  ^  2.  It  is  reasonable  to  assume  in  this  case  that 
n 
r 

the  observers  lie  within  the  relaxation  zone.  Making  this  zone  infinite 
in  size,  one  finds  that  [u]N  =  0  (1/w)  for  w  <<  1.  This  merely  expresses 
the  fact  that  a  net  static  displacement  is  to  be  observed. 

-3 

Both  of  the  extreme  stress  conditions  yield  average  slopes  of  w  for 
the  high  frequency  end  of  the  spectrum.  This  fact  is  critical  when  one 
tries  to  explain  data. 

A  major  effort  has  been  undertaken  to  add  these  new  theoretical 
developments  to  existing  computer  programs.  Numerical  computations  constitute 
the  most  efficient  method  of  determining  the  effects  of  the  different  source 
parameters  on  the  radiation  field  and  the  asymptotic  behavior  described 
above  provides  an  important  check  on  the  accuracy  of  the  numerical  codes. 


m!Vl  yj*urwi  w  WWW  V/j  V"r^“'  jju-irww^! 
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(2)  Surface  Waves  from  a  General  Elastodynamic  Source 
In  a  Vertically  Inhomogeneous  Half-space 

As  our  source  in  a  locally  homogeneous  region,  we  take  the  slightly  modified 
elastodynamic  source  form  of  Archambeau  (1968). 
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where  <J>g  and  if;  are  the  Fourier-time  transformed  compressional  and 
Cartesian  shear  potentials  (j=l,2  and  3)  respectively.  In  order  to  express 
these  potentials  in  terms  of  the  separable  solutions  to  the  Helmholtz 
equation  in  cylindrical  coordinates,  we  use  the  following  relation 
(Harkrider  and  Archambeau,  1973) 
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-iv  | z— h | 

e  V  Jm(kr)  kdk 


where 


V  =  kr 
v  v 


(k2  -  k2>1/2 

V 


-i(k2-k2)i/2 


;  k  <  kv  ,  kv  -  to/ v 
;  k  >  k 

v 


p“(D  s  d-Om/2  P„m)(V 

Tj(0  =  (?2-i)m/2  P^m)  (O 

P(m)(?)  *  —  P  (O 
n  w  ,_m  n 
dt, 


v  is  either  a  or  8,  the  compression  or  shear  velocity  respectively  and 
(r,z)  =  (0,h)  is  the  source  location. 


Making  use  of  this  relation,  we  can  rewrite  equations  (1)  as 


sin  m4> 


-lkr  z-h 

a1  1 
e 


r 

a 


J  (kr)  dk 

m 


(3) 


Next  we  obtain  an  expression  for  the  cylindrical  SV  potential, 
which  is  a  convenient  potential  for  our  cylindrical  coordinate  system,  in 
terms  of  the  Cartesian  SV  potentials  given  in  equation  (3).  The  vertical 
displacement  integrand,  w,  of  its  k  integral  is  related  to  the  compressional 
and  Cartesian  SV  potential  integrands  by 
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.ii  *  .  !!i 

W  3z  3x  3y 


and  in  terms  of  the  congressional  and  SV  potential  integrands  by 


.  -  It  +  iA 


which  by  inspection  yields  the  relation 


i_  (0*2  _ 

2  \  3x  3y  / 


Performing  the  above  operation  and  comparing  with  the  cylindrical  SV 
potential 


-§/{ 


E  cos  m(j)  +  F  sin  md>  >  - 

m  m  (  r. 


-ikrg|z-h| 


J  (kr)  dk 
m 


we  obtain  the  following  relation  between  coefficients  as  derived  in 
(Harkrider  and  Archambeau,  1973). 

2kIm  -(»^1  +  »i“) 


where 


C(j)  and  are  zero  for  m  >  n  and  m  <  0  , 

m  m  ’ 
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and  in  addition 


C<2)  =  D<X> 


c(1)  =  -  n(2) 

C0  D0 


Fq  =  0. 


The  cylindrical  source  potentials  given  by  equations  (3)  and  (8)  may 
now  be  substituted  into  the  multilayer  formulation  of  Harkrider  (1964).  But 
first  we  note  that  alternating  terms  in  the  infinite  series  in  equations  (4) 

are  of  opposite  sign  depending  on  whether  z  is  greater  or  lesser  than  h. 
We  separate  the  series  such  that 


Am  =  A®  +  A° 
m  m  m 


where  the  e  superscript  denotes  a  new  series  made  up  of  the  terms  with 
m+n  even  and  the  0,  a  series  formed  by  terms  with  m+n  odd.  A  similar 
separation  is  done  for  the  other  source  coefficients.  The  new  coefficients 
have  the  following  property 


AJz  >  h) 


Am(z  <  h> 


A“(z  > 
m 


-  A^(z  <  h) 


Defining 


40. 
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c(3)o 
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i2k 


5(3)° 

m 


-  i2k2y  C(3)e 
in 


-  i2k2y 


D(3)e 


Following  Harkrider  (1964)  we  obtain  as  our  integral  solution  for  the 
vertical  displacement  at  the  surface  of  our  inhomogeneous  half-space 


<Wo  > 


^  j  ^1  ^RlltA]m  ~  R12[B]m  +  R13Zm 


J  (kr)  dk 
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where 


Fe  -  -  Ru  -  [t]R12 


[T]  = 
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[A]m  '  [-'V43  \  +  <V33  K  -  <V23  {Z„  +  <Vl3  {X„] 

[B,m  =  [(Ars)42  ®m  '  (\s' 32  ™m  +  (ARS>22  5Zn,  '  <ARS,12  6Xm  1 

J  (15) 

and 

Zm  -[-Vu  5Um  +  <V31  4Hm  -  (ARS>21  “m  +  <Wll  «.] 

(For  symbols  used  without  definition  here,  refer  to  Harkrider  (1964),  (1970) 
and  Harkrider  and  Flinn  (1970) . 

The  matrix  as  defined  in  Harkrider  (1964)  is  the  layer  product 
matrix  which  gives  the  displacement-stress  vector  associated  with  P-SV  motion 
at  source  depth  in  terms  of  the  surface  displacement-stress  vector.  The 
integral  solution  given  by  equation  (14)  is  also  valid  for  a  vertically 
inhomogeneous  half-space  where  is  the  linear  transformation  of  the  P-SV 
displacement  stress  vector  from  the  surface  down  to  the  source.  The  only 
restriction  on  this  form  of  the  solution  is  that  at  some  depth  the  media 
is  terminated  by  a  homogeneous  half-space. 

The  surface  azimuthal  displacement  due  to  SH  waves  is  given  by 


<Dm(2) 
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dJ  (kr) 
in 


d(k  ) 
ra 


dk 


(16) 


fl  =  -  (V21  ■  (AL>n^  rf 


u 


N(1)  =  i 


m 


_(AL}22  "  (Vl2yi  rgJl][(ALS)ll  6Ym  “  (ALS)12  5VmJ 


(17) 
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The  \s  113  the  transfer  matrix  for  the  displacement  stress  vector  associated 
with  SH  motion  down  from  the  surface  to  the  source  depth.  is  the  transfer 
matrix  from  the  surface  down  to  the  depth  at  which  the  terminating 
homogeneous  half-space  begins  with  elastic  properties  denoted  by  subscript  l. 

Evaluating  the  residue  contributions  of  equations  (14)  and  (16) ,  in 
order  to  obtain  the  surface  displacements  due  to  Rayleigh  and  Love  waves 
respectively,  yields 


and 


■>.  Ajhiu-etear.  slid  i-jia  ■ ..  ■-  r .  ■.■in  •  ■  u  l-i-  L.  'A  .•  t.  o-i. 
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Using  equations  (13)  the  solutions  can  be  written  as 
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E  (^3)e  cos 

m=0  \ 


+  D(3)e  sin 
m 


—  (kr) 


and  the  o  superscripted  variables  defined  similarly. 


(3)  Hear  Field  Wave  Propagation  in  Layered  Media  using  SlmnlP 
Dislocation  and  Explosion  Sources . 
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We  begin  by  first  reviewing  the  treatment  of  a  spherically  symmetrical 
source  followed  by  the  general  treatment  of  dislocation  models.  Our 
approach  is  to  develop  the  mathematical  tools  for  a  very  simple  model, 
an  infinite  homogeneous  medium,  that  can  be  used  in  treating  the  "n"  layer 
model  after  the  application  of  generalized  ray  theory.  To  apply  generalized 
ray  theory  we  must  reduce  all  sources  to  simple  displacement  potentials 
of  the  P,  SV,  and  SH  type.  This  is  easily  accomplished  for  symmetric 
sources  like  explosions  but  difficult  for  earthquake  sources.  We  will 
present  such  potentials  later. 

The  radial  displacement  for  a  simple  point  source  can  be  expressed  in 

terms  of  the  potential 

a)  uR(R,t)  -Ij+ct  -f> 

where 

(2)  <j)(R,t)  =  4jq  6(t  -  R/V)/R. 

R  is  the  radial  coordinate  and  V  is  the  velocity.  The  parameter  <|>  is  a 

o 

constant  with  units  of  volume  (assumed  to  be  unity)  and  we  are  assuming  a 

delta  function  time  history.  Taking  the  Laplace  transform  of  (2)  yields 

R 

_  "  V  8 

(3)  0»(r,s)=£__  . 

Applying  the  Sommerfeld  transformation  we  can  write  J  in  terms  of  cylindrical 
coordinates 
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where  k  is  the  horizontal  wave  number.  Next  we  go  from  J  to  K  (modified 

o  o 

Bessel  function)  using  the  basic  Cagniard-de  Hoop  notation 


k  =  -  isp 


and 


(5) 


where 


-s  nJz| 

KQ(spr)  e  dp 


Now  we  must  take  the  inverse  Laplace  transform,  in  general  we  will 
be  working  with  K^spr)  so  we  will  allow  n  to  be  unspecified  for  the 
moment.  Equation  (5)  can  be  solved  exactly  but  there  are  some  useful 
approximations  to  discuss.  Using  the  asymptotic  expansion  (9.7.2  Abramowitz 
and  Stegun)  we  have 


(6) 


Kn(x) 


1  + 


(y-1) (y-9)  . 
2(8x) 2  + 


where  p  =»  4n2.  The  series  converges  the  fastest  for  small  values  of  n. 
However,  assuming  n  =  2  (the  highest  order  needed  in  dislocation  modeling) 
and  x  *  3  we  get 

M3)  -  e'3  £  1  +  .623  +  .091  j  . 

This  approximation  is  98 7.  accurate. 


'  - - 


■  . . 
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If  we  assume  large  (spr)  we  have 
/  \l/2 

(7)  Ko(Spr)^^j  e'spr 

and  substituting  into  (5)  we  obtain 


. . 

which  can  be  treated  by  line  source  theory. 
Following  the  de  Hoop  transformation 


dp 


t  =  pr  +  rys  z  >  0 

we  obtain 


The  asymptotic  solution  becomes 


> 


l 
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Equation  (10)  is  to  be  evaluated  along  the  contour  defined  by  (9).  A 
further  simplication  of  (10)  can  be  obtained  by  making  a  so  called 

_  J  _  L 

"first-motion"  approximation.  For  values  of  t  near  R/V,  p  «  — — — 
and  nv  ■  COy-  —  and  (10)  reduces  to 


(13)  4>  (r,z,t)  «  I  &(t  -  , 

which  was  the  starting  equation  (2).  Note  that  p  »  9iy  h  (the  ray 
parameter)  in  the  first-motion  approximation.  This  approximation  is 
useful  for  comparing  Cagniard  solutions  of  dislocation  models  with 


conventional  far-field  results. 
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The  exact  solution  of  (5)  can  be  obtained  by  applying  transform 
theory.  We  note  that 


and  applying  the  shift  rule 


50. 


The  integral  for  various  values  of  t  can  be  evaluated  using  the  transfor¬ 
mations  used  by  Helmberger  (1968). 


Earthquake  sources 


In  this  section  we  will  examine  some  relatively  simple  dislocation 
models.  Starting  with  Haskell's  representation  for  shear  faulting  it  is 
relatively  easy  to  devise  the  displacements  for  double  couples  in  an 
infinite  medium.  The  solution  for  a  strike-slip  fault  becomes 


77,  v  „  T  3  /  92A  1  3A\ 

w<r-z’“>  -K  [i?  (1^-7  a) 


2^/92a  _  1  9A  \ 
r  \  3r2  r  3r  J 


sin  26 


V(r,z,oj)  =  K 


+  ks  3?  cos  26 


-Krf(|4-I|i)+  k2|B  ] 

n  ’  ’  3r  \  3r2  r  3r  /  3  3r 


sin  26 


where  W,  V  and  q  are  the  displacements  in  the  vertical,  azimuthal  and  radial 

coordinatee.  The  parameters  are 

.  R  .  R 
“1“  T  “iw  ~ 

_  p  „  a 

.  e  e 


B  *  e 


-i»i 
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K  -  -  h  LKp(“)  *  m  a 

4ttu)2p  ’  Kg  3 


and 


D(w)  is  the  Fourier  transform  of  the  time  history  across  the  fault. 

We  are  assuming  an  instantaneous  motion  over  a  rectangle  with  dimensions 
L.H,  It  is  easy  to  generalize  to  finite  rupture  velocity  by  adding  the 
moving  source  directivity  (see  Mikumo,  1972).  We  will  add  such  features 
later  along  with  multiple  sources  as  well.  Our  approach  is  to  compute  and 
understand  the  three  basic  faults,  strike-slip,  dip-slip  and  45°  dip 
slip.  Burridge  and  Knopoff  (1964)  show  that  a  linear  combination  of  these 
three  solutions  will  represent  any  fault  orientation.  We  will  treat  the 
strike-slip  dislocation  as  an  example  although  any  motion  can  be 
handled  using  our  approach. 

Next,  we  breakdown  these  displacements  in  terms  of  potentials, 

_  =  ii  +  l5+k2  r 

W  3z  3z2  6  V 


T7  =  I  +  1  d  'P  _  lA 

r  30  r  3z30  3r 


3r  3r  3z  r  38 


After  some  effort  the  potentials  become 


-  K  /  k2Fa  J2(kr)  sin  29  dk 


—  /  Fg  J2(kr)  sin  20  dk 


“if 


Fg  J2(kr)  cos  20  dk 
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where 


-<k^)1/2|z| 


k  e 


(19) 


F„  = 


‘V  -'.2  v2n1/2 


(k2-  K2)' 

We  are  concerned  with  the  evaluation  of  the  integrals  in  (18),  so 
we  need  only  consider  the  field  variable: 


(20)  C(r,z,w)  = — /  Fy  J2(kr)  dk 


-f 


Changing  variables 


u)  =  -  is  and  k  =  -  isp 


we  obtain 


(21)  C(r, z, s)  =  |  s  Im  j  K2(spr)  e  '  dp 
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i°° 


•SHyZ 


0 


where 


nv  =  (  W  '  p2 


1 1/2 


Equation  (21)  is  almost  identical  to  (5)  of  the  previous  section  and 
can  be  solved  following  the  same  technique.  Since  K2  decays  slower  than 
Kq  we  know  that  q  will  develop  a  tail.  That  is,  the  polarity  of  the  second 
term  of  the  series  given  in  equation  (6)  is  negative  for  n  =  0  and 
positive  for  n  =  2  indicating  a  long  period  enhancement.  Suppose  we 
express  c,  in  a  series 

C  =  Ci  +  C2  +  ?3 


where  the  £'s  indicate  the  various  terms  in  this  expansion; 


This  series  can  be  readily  evaluated  and  is  applicable  for  all  periods 
such  that 


T  < 


2-rrr 

3V 


For  small  distances  and  long  periods  we  must  use  expressions  similar  to 
(16).  Before  generalizing  to  a  layered  half  space  we  will  consider  the 
first-motion  approximations  of  the  equations  given  by  (18).  First  wo  note 

^  ~  /"1  =  R<5^t  V } 


that 


where 
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sinh 

V 


COST 


and  nv - y — ‘  >  as  indicated  from  the  first  section.  We 


can  then  write 


*(r,2,s)  --  +  W4^P*;>  (32p2)  L 
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(23) 
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o _ 0_ 

4tt 


(^)2 


H(t-R/B) 

R 


sin  2 


assuming  D(s)  -  Dq/s  ,  that  is  we  assume  a  sti>p  function  time  history  of 
the  displacement  across  the  fault.  The  moment  is  just  ji(LHD  ). 


Similarly, 


Bo(LHno>  h  \  H(t- 


R/V) 


(24)  x(r,2,t)  -  -  "  -(fr) C0S  20 


”1  s  — 

It  is  convenient  to  define  a  new  potential  =  —  ip  which  reduces  to 


/ocx  ,1/  n  Bo(LHDo}  /  cos  h  \  H(t-; 

(25)  ^(r.z.t)  =  -  — -  ^  -jr-  )  — 


-R/B) 


sin  20 


where  we  used  n„  = 


cos  h 


3  B 

These  first-motion  approximations  can  be  compared  with  the  results 
of  geometrical  ray  theory  by  computing  the  displacement  due  to  P  and  SV  and 
SH  waves  where 


UP  •  +  i*l 1 


(26) 


(27) 


“sv  '  (“,2  +  $ 


U  =  V 
SH  x 


For  example,  note  that  the  operations  can  be  written 


L  1  c(r,z,s)  =  -  L  1  ^(spK  j 


. -i  r  a  ,  .i  ,-ir  i 
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and  using  the  definitions  (17)  we  have 


sin  20 


-(g)(5) 


R(h ,9)  5 (t-R/a) 
R 


where 


R(h,0)  = 


sin2h  sin  28 


This  is  the  far-field  expression  for  the  field  given  by  iJen-Menahem, 

Smith  and  Teng  (1965). 

The  form  of  the  solution  to  use  is  dictated  by  the  value  of  r,  the 
source  time  history  and  the  model  response.  It  would  appear  that  all  three 
solutions,  exact  integral  evaluation  given  by  (16),  the  power  series 
development  (22),  and  the  first-motion  approximation,  have  domains  of 
interest.  We  are  primarily  interested  in  the  first  two,  especially  the 
power  series  solution  for  application  to  the  near  field  radiation  field 


observations . 


56. 


Generalization  to  layered  half  space 

Using  the  method  of  generalized  rays  we  can  construct  solutions 
by  the  same  technique  as  used  for  the  case  of  no  azimuthal  dependence. 
We  will  not  list  all  the  expressions  here  but  we  will  work  out  the 
equations  for  the.  ^-potential  as  an  example.  For  an  incident  P-wave 
at  the  free  surface  we  have 


7m  J p2  Fa(t)  G^(z)  sin  20  dp 

where 


and  the  time  history  across  the  fault  is  assumed  to  be  a  step  function 
in  displacement.  The  sources  is  situated  at  zq.  There  will  also  be  a 
converted  SV  wave  which  can  be  written 
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If  one  substitutes  these  equations  into  the  stress  equation,  the  zero 
stress  conditions  at  the  boundary  are  satisfied.  The  vertical  displacement 


becomes 
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w  (r,z,e,s)  =  Mq  -  s  Im  /  p2  M(ct )  RPZ  sin  20  dp 


where 
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K1 <sPr)  +  —  K2 (sPr) 
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2r,c  (Pb  “  P2) 
62  R(p) 


R(p)  -  (n2  -  p2)  +4  p2n  n 

3  a  t 


R(P)  is  just  the  RayleiSh  denominator.  The  radial  displacement  produced 


uy  the  <p  potential  becomes 


Mo(r)s  uf 


(P2)  M(a)  RPR  sin  20  dp 


where 


RPR  = 


4  P"c"b 
e2R(p) 


There  will  also  be  a  tangential  component  generated  by  this  tj>  potential 
since  it  contains  a  "0"  dependence, 


'2-e's:>  ■  Mo(f)  r  J 


(p2)  RPT  (2  cos  20)  M  (r)  dp 


Rpx  -  1  +  Rpp  -  BnBRre. 


where 
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>rhe  reflection  coefficients  Rpp  and  Rps  are  given  in  Helmberger  (1968). 
The  various  receiver  functions  RPZ ,  RPR  and  RPT  all  contain  the  Rayleigh 
denominator  and  thus  there  will  be  a  Rayleigh  wave  propagating  along 
with  motion  in  all  three  coordinates.  These  solutions  can  be  transformed 
back  to  the  time  domain  by  breaking  the  solution  into  far  and  near  fields 
denoted  Wf-^  and  WN^  ,etc.  Applying  the  theory  developed  in  previous 
sections  we  obtain 


WF 
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and  for  the  near  field  we  obtain 
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It  is  useful  to  examine  the  high  frequency  solution  by  expanding  M(V) 
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and 
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The  second  term  is  simply 
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The  interpretation  of  these  equations  is  relatively  simple.  If 
one  neglects  p2  RPZ  in  JW^1  ,  we  obtain  for  W1^  a  delta  function  time 
response  multiplied  by  the  source  strength  M  and  divided  by  the 
distance  traveled  .  The  time  function  is  the  derivative  of  the  source 
time  function.  The  p2  factor  is  a  source  correction  for  take-off  angle  as 

mentioned  earlier.  The  RPZ  function  is  equal  to  2  at  vertical  incidence, 
which  is  the  plane  wave  surface  effect  (Knopoff,  Fredericks,  Gangi  and 
Porter  (1957)).  This  function  yields  the  Rayleigh  wave  at  p  =  — ,  where 
VR  is  the  Rayleigh  velocity.  The  time  function  will  look  like  the 
derivative  of  a  Gaussian  function.  The  second  term  W?  will  also 
contain  a  Rayleigh  wave  but  will  appear  as  a  Gaussian  and  decay 
l / r  faster. 


TH-  v-V! 
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The  displacements  produced  by  the  other  potentials  x  and  i p  can  be 
computed  following  the  same  procedure.  A  complete  presentation  in  matrix 
form  will  be  given  later.  To  generalize  to  "n"  layers  we  need  only  apply 
the  concepts  of  generalized  reflections  and  transmission  coefficients. 
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IV.  Summary  of  Theoretical  Source  Predictions:  Implications  for 
Discrimination  of  Earthquakes  and  Explosions 

In  this  section  the  preliminary  applications  of  the  theoretical  results 
of  section  III  are  summarized  in  terms  of  their  relevance  to  discrimination 
of  earthquakes  and  explosions.  Additional  work  under  a  separate  program 
(AFOSR  Contract  F44620-72-C-0078,  Seismic  Phenomena  Connected  with  Earthquakes 
and  Explosions)  aimed  at  a  systematic  modeling  of  both  earthquakes  and 
explosions  for  purposes  of  defining  and  understanding  discrimination  criteria 
is  also  underway  and  to  some  extent  the  work  reported  in  this  section 
overlaps  work  under  this  second  program.  However,  the  work  here  is  distinct 
in  the  sense  that  we  have  concentrated  on  the  near  field  and  the  more  complex 
problems  related  to  observations  in  this  range  while  in  the  former  program 
we  are  mainly  concerned  with  the  far  field  "teleselsmlc"  radiation.  Nevertheless 
our  investigations  of  the  near  field  and  the  attempt  to  incorporate  complex 
source  models  within  the  various  wave  propagation  programs  for  studies  in  the 
near  field  range  have  resulted  in  a  more  complete  understanding  of 
earthquake  radiation  and  the  underlying  basis  for  discrimination. 


A 
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In  particular,  theoretical  methods  were  devised  which  enabled  us  to  express  the 
complicated  source  models  of  Section  IH-1.  above  as  equivalent  point  (multipolar) 
sources.  (These  methods  apply  to  both  analytical  source  models  and  numerical 
source  models.)  The  approach  was  first  utilized  to  combine  explosive  source 
models  and  the  tectonic  release  effects  together,  and  then  to  incorporate 
this  composite  source  into  a  wave  propagation  program  for  both  surface  and 
body  waves  (ray  theoretic  approximation)  so  that  we  are  able  to  predict  tee 
teleseismic  radiation  of  first  and  later  arrival  body  waves  (including  surface 
reflections)  and  Love  and  Rayleigh  surface  waves  at  any  point  in  or  on  the 
earth.  The  same  procedure  was  followed  for  earthquake  sources  in  the  present 
program,  wherein  a  variety  of  complex  relaxation  source  models  arc  included 
and  for  which  the  teleseismic  radiation  fields  of  body  waves  and  surface 
waves  can  be  predicted  at  any  spatial  point  and  in  either  the  frequency  or 
time  domains.  The  very  near  field  wave  propagation  is  handl.  d  using  the 
Cagniard  method.  This  approach  gives  results  directly  in  the  time  domain. 


Spectra  can  also  be  obtained  from  these  results. 

Using  this  predictive  capability,  we  began  a  systematic  investigation 
of  the  radiation  fields  from  the  variety  of  source  representations 
described  in  Section  III  in  order  to  determine  the  spectral 
characteristics  of  tcleseismic  signals  from  these  sources  and  to  compare 
these  spectral  predictions  with  observations.  Our  purpose  was  to  obtain  <- 
best  model  (or  model  series)  for  earthquakes  using  data  from  events  in  vr.r 
environments.  Secondly,  we  wanted  to  determine  the  spectral  details 
number  of  particular  seismic  phases,  in  particular,  the  first  arriving 
pP  phases  (separately  and  combined)  and  the  Love  and  Rayleigh  surface  w.-.ve 
for  these  "best"  models.  Some  of  the  more  important  results  were  as  fo  " 
(a)  Earthquakes  appear  to  be  reasonably  well  modeled  by  relaxation 
source  models  in  that,  of  the  rather  incomplete  set  of  spectral 
observations  made,  we  can  obtain  a  reasonable  fit  to  any  of  the 
data.  In  particular,  either  flat  appearing  or  strongly  peaked 
appearing  spectra  can  be  fit  by  the  model  by  choice  of  a  larger 
prestrained  region  or  small  prestrained  region  relative  to  the 
failure  or  fault  zone  dimension,  so  that  quite  different  looking 
spectral  observations  can  be  fit  by  the  same  basic  model.  Indeed, 
the  model  actually  predicts  that  a  rather  wide  range  of  possibilities 
for  the  radiation  field  spectra  are  to  be  reasonably  expected.  How¬ 
ever,  we  will  have  to  wait  for  additional  data  to  be  sure  that  this 
range  is  actually  met  in  nature;  we  will  conclude  that  earthquakes 
within  this  predicted  range  are  possible  only  if  we  are  sure  the 
model  is  in  close  agreement  with  a  wide  range  of  data. 


(b)  The  P  and  S  wave  displacement  spectra  predicted  by  relaxation 
sources  and  partially  verified  by  observations,  are  such  that  the  P 
and  S  wave  spectra  have  characteristic  frequencies  fP  and  fS  which 

c  c 

effectively  divide  the  spectra  into  high  and  low  frequency  parts. 

The  characteristic  (i.e.,  comer  or  peak)  frequencies  are  soim;i 

azimuth  dependent  and  f®  <  fP  generally,  although  this  relation:.;.  ■ : . 

can  be  reversed  in  some  azimuths.  The  behavior  at  high  frequency 
P  s 

^  >  ^or  suc^  that  the  average  amplitude  of  the  spectrum 

decreases  rapidly  from  a  maximum  level  at  fP(f^)  reaching  a  mean 
slope  proportional  to  1/f3  for  f  »  f£(f®).  Superposed  on  this 
average  behavior  are  secondary  maxima  and  minima  that  are  due  to 
rupture  propagation  effects  resulting  in  constructive  and  destructive 
inference.  At  low  frequencies  f  <  f^(f^)  the  spectrum  is  even  more 

complex  in  that  both  near  and  far  field  effects  are  present,  and 
spectral  shapes  at  long  periods  depend  on  the  distance  of  the  obser¬ 
vation  point  from  the  source.  The  near  field  effects  are  dominated 
by  terras  which  behave  as  1/f  at  low  frequency  and  this  term  can 
dominate  within  the  entire  range  f  <  fP(f^)  at  near  source  distances. 
At  larger  distances , the  far  field  terms,  having  different  distance 
dependence,  dominate  and  the  spectral  behavior  of  this  far  field 
component  is  dependent  on  the  dimensions  of  the  prestress  region. 

In  particular,  the  spectrum  is  flat  to  zero  frequency  if  the  pre¬ 
stress  region  is  taken  to  be  infinite  (which  is  impossible  in  the 
earth,  of  course),  and  decreases  for  f  «  fP(f^)  for  the  prestress 
region  finite  with  a  slope  of  f2  for  f  small.  For  a  small 


prestress 


zone  (prestressed  zone  dimension  comparable  to  the  rupture  dimension) 
the  spectrum  is  strongly  peaked  at  f£(fjj)  decreasing  immediately 
toward  the  lower  frequency  end  of  the  spectrum  with  the  f2  slope 
while  for  a  larger  prestress  zone  dimension,  the  peak  broadens  and 
as  the  prestress  zone  increases  the  peak  continuously  broadens 
being  essentially  flat  until  a  frequency  near  f'R(f'S)  where  it 
begins  to  decrease  again  with  the  f2  slope.  In  general,  the  S  wave 
spectra  (either  SV  or  SH)  are  from  5  to  10  times  larger  in  magnitude 
than  the  P  wave  spectra  and  shifted  somewhat  to  lower  frequency  but 
otherwise  of  similar  overall  shape. 

(c)  The  spectral  characteristics  discussed  in  (b)  above  have  been 
found,  theoretically,  to  scale  in  a  simple  way.  That  is, a  complete 
i3L  spectrum  for  an  earthquake  model,  with  failure  parameters 

specified  in  terms  of  prestress  magnitude  and  orientation  a(0)  f nuj  t 

ij  ’ 

or  failure  zone  dimension  L,  rupture  velocity  VR,  intrinsic  material 
velocities  Vp  and  Vg ,  and  relaxation  zone  or  initial  prestress  zone- 
characteristic  dimension"  R^  ;  can  be  scaled  to  yield  the  spool,  r.w 
for  any  sized  earthquake  of  similar  type.  These  scaling  laws  _ 

to  be  in  at  least  rough  agreement  with  observations.  The  scaling 
law  is  as  follows: 

(1)  The  characteristic  frequency  for  P  waves: 

fc“  W/1 

p 

In  fact,  fj,  is  essentially  equal  to  VR/L  except  that  it  varies 
somewhat  with  azimuth. 


The  characteristic  frequency  for  S  waves: 


4“  (vv  ^v'1 

s 

Again  f^  is  very  nearly  equal  to  this  quantity  except  for 
azimuthal  variations  which  are  somewhat  larger  than  those  for 
the  P  wave  radiation  (e.g.  f  '  can  be  larger  than  f^  at  some 
azimuths)  . 

(2)  The  spectral  amplitude  | u^ |  scales  with  rupture  dimension 
L  as : 

|u  |  a 
P 

Similarly  for  S  waves  : 

| u  |  a 
s 

The  shape  of  the  spectrum  is  not  altered  with  changes  in  L, 
holding  all  other  parameters  fixed.  (Note,  however,  that  as 
L  changes  both  fc  and  shift.) 

(3)  The  spectral  amplitudes  scale  directly  with  the  prestress 
(or  prestrain,  e^^  more  properly),  so: 

i%l  a  le^l  and  \%\  «  |e^}| 

(4)  The  width  of  the  spectral  peak  (or  flat  portion  of  the 
spectrum  which  begins  at  f£  (or  f^  for  S  waves)  is  controlled 

*  p  ip 

by  Rg.  The  frequency  f^,  (f^,  )  at  which  the  spectrum  falls  off 

2 

with  decreasing  f  (as  f  )  scales  as: 
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Therefore  the  spectral  width  of  the  peak  (the  flat  portion  of  ..nc 

spectrum  is  Af  =  f^  -  f^  or  Af  =  f^  -  f^°.  Again  the  frequencies 
1 P  1  S 

an<^  aPPear  P°  be  nearly  equal  to  the  expressions  given 
above  (the  factor  of  proportionality  is  nearly  unity  but  varies 
with  azimuth  somewhat) . 

(d)  The  radiation  patterns  for  P  and  S  waves  correspond  to  superposed 
multipole  patterns,  but  the  quadrupole  is  dominant  at  all  frequencies 
(double  couple  force  equivalent) .  However,  at  high  frequencies  the 
higher  order  multipoles  became  significant  and  the  resulting  observed 
patterns  are  distorted  quadrupole  (4  lobe)  patterns,  usually  with 
high  amplitudes  in  the  direction  of  rupture  propagation.  These  high 
frequency  effects  are  due  to  rupture  propagation  (moving  source 

effects) .  The  patterns  are  nearly  pure  quadrupole  at  frequencies 

P  S 

less  than  f  and  f  . 

I-#  c 

In  view  of  the  previously  summarized  results,  we  considered  the 
implications  of  these  source  properties  insofar  as  discrimination  was 
concerned.  The  discriminates  considered  in  some  detail  so  far  were  nq  vs. 

Mg,  and  P  wave  spectral  shapes.  The  -  M  results  for  earthquake  can  be 
summarized  as  follows : 

(a)  and  Mg  were  estimated  using  the  log  of  the  P  wave  amplitude  at 
1  cps  for  and  the  log  of  the  SV  amplitude  at  .05  cps  for  Mg.  One 
standard  event  was  computed  and  the  scaling  laws  previously  described 


were  used  to  obtain  results  for  larger  and  smalJer  events.  (This 

will  be  replaced  later  by  actual  time  domain  measurements  using 

synthetic  seismograms  generated  from  the  spectra,  and  employing 

tne  actual  field  procedure  followed  in  obtaining  m^  and  M  ).  for 

earthquakes  the  m^  vs .  curve  had  a  one  to  one  slope  at  low 

magnitudes.  Because  of  the  expected  variations  in  prestress  a.  '! 

dimensions  of  the  zone  ci  prestress  we  actually  obtain  a  band  i:i 

the  m^  -  plane  within  which  earthquakes  should  occur.  The  b.rr.u 

width  is  about  one  unit  on  the  (vertical)  M  scale.  The  band  uhen 

s 

has  the  45  degree  slope  in  the  -  Ms  plane.  Further  at 
magnitudes  which  appear  to  correspond  to  from  5+  to  6.5+,  the  curve 
quite  abruptly  becomes  vertical  indicating  a  maximum  rr^  of  around 
6.5.  Larger  fault  lengths  involving  larger  earthquakes  will  not 
give  m^  any  larger  than  about  6.5.  But  witn  very  high  prestress 
of  the  order  of  kilobars ,  this  maximum  m^  could  be  as  high  as  7.  The 

curve  (or  band  more  properly)  continues  upward  until  at  around  M  of 

s 

11)  it  terminates.  (No  earthquake,  however  great  in  length,  would 
have  a  larger  M  than  11,  for  fl  measured  in  this  way.)  The  fault 
lengths  corresponding  to  the  critical  points  for  earthquake  vs.  M 
are  near  10  km  (for  the  cutoff)  and  near  900  km  (for  the  Ms  and 
cutoff  or  termination  point).  The  earthquake  population  should  lie  in 
c.  band  of  one  magnitude  unit  (roughly)  and  this  arises  from  the  possibl 
variations  among  earthquakes  of  the  rupture  velocity,  the  prestress, 
and  dimension  of  the  prestress  zone. 


